        subroutine ROHF(NORB,X_T,F_a,F_b,P_a,P_b,S,Hcore,U,n_el_a,n_el_b,E_nu)
        !!!F_a F_b P_a P_b is convergent in UHF;S is overlap matrix
        implicit none
        integer :: i,j,k,l,m,count_iter
        integer :: NORB,n_el_a,n_el_b
        real*8  :: U(NORB,NORB,NORB,NORB),Hcore(NORB,NORB)
        real*8  :: X_T(NORB,NORB),X_A(NORB,NORB)
        real*8  :: F_a(NORB,NORB),F_b(NORB,NORB)
        real*8  :: F_RO(NORB,NORB),F_ROHF1(NORB,NORB)
        real*8  :: F_ab(NORB,NORB),eye(NORB,NORB)
        real*8  :: C_RO(NORB,NORB),S(NORB,NORB)
        real*8  :: P_v(NORB,NORB),P_t(NORB,NORB)
        real*8  :: P_a(NORB,NORB),P_b(NORB,NORB)
        real*8  :: P_a1(NORB,NORB),P_b1(NORB,NORB)
        real*8  :: pv(NORB,NORB),pc(NORB,NORB)
        real*8  :: temp(NORB,NORB),tem(NORB,NORB)
        real*8  :: eig(NORB,NORB),eig_vec(NORB,NORB)
        real*8  :: orb_e(NORB)
        real*8  :: E_ele,E_ele1,delta,delta1
        real*8  :: E_nu
        logical :: logical_e
        !!!inital guess
        write(21,*)'****ROHF program activate****'
        if (n_el_a>n_el_b)then  !!! Calculate the O2 triplet;n_el_a=9;n_el_b=7
                write(21,*)'ROHF Correct distribution'
        else
                write(*,*) 'ROHF Error distribution '
                stop
        endif
        !!!*************************************************************
        !!Construct ROHF Fock matrix
        !!ref. J. Chem. Phys. 133, 141102 (2010)!!
        !======== ======== ====== =========
        !space     closed   open   virtual
        !======== ======== ====== =========
        !closed      Fc      Fb     Fc
        !open        Fb      Fc     Fa
        !virtual     Fc      Fa     Fc
        !======== ======== ====== =========
        !write(*,*)S
        F_RO=0.0d0
        eye=0.0d0   !!!Construct E matrix
        do i=1,NORB
                eye(i,i)=1.0d0
        enddo
        
        !!!Construct F_RO matrix
        P_t=P_a+P_b             !!! D^a + D^b
        P_v=P_a-P_b             !!! D^a - D^b
        pc =matmul(S,P_t)       !!! S*(D^a + D^b)
        pv =matmul(P_v,S)       !!! (D^a - D^b)*S
        F_ab=F_a-F_b            !!! F^a-F^b
        tem=eye-pc              !!! (1-S*(D^a+D^b))
        temp=matmul(matmul(tem,F_ab),pv)          !!! temp=(1-S(D^a+D^b))(F^a-F^b)(D^a-D^b)S
        F_RO=0.5d0*(F_a+F_b+temp+transpose(temp)) !!! F^RO=1/2*(F^a + F^b + temp + temp^T)
        !!! End Construct!!!!

        X_A=transpose(X_T) 
        F_ROHF1=matmul(matmul(X_A,F_RO),X_T)      !!! FC=SCe--> F'C'=C'e 
        eig_vec=F_ROHF1
       
        call DIAG(eig_vec,eig,NORB)          !!! eig_vec matrix output eigenvector 
        C_RO=matmul(X_T,eig_vec)             !!! C=XC'    
        call Density1(C_RO,P_a,NORB,n_el_a)  !!! inital guess D^a
        call Density1(C_RO,P_b,NORB,n_el_b)  !!! inital guess D^b 
        call E_ROHF(NORB,F_a,F_b,P_a,P_b,Hcore,E_ele) !!! obtain E_ele
        E_ele1=E_ele
        logical_e=.true.
        count_iter=0

        !write(*,*)'111'
        do while(logical_e)
        !do m=1,100
        count_iter=count_iter+1
        F_a=Hcore;F_b=Hcore
        do i=1,NORB
          do j=1,NORB
             do k=1,NORB
                do l=1,NORB
                   F_a(i,j)=F_a(i,j)+(P_a(k,l)+P_b(k,l))*U(i,j,k,l)-P_a(k,l)*U(i,l,k,j)
                   F_b(i,j)=F_b(i,j)+(P_a(k,l)+P_b(k,l))*U(i,j,k,l)-P_b(k,l)*U(i,l,k,j)
                enddo
             enddo
          enddo
        enddo
        P_t=P_a+P_b             !!! D^a + D^b
        P_v=P_a-P_b             !!! D^a - D^b
        pc=matmul(S,P_t)        !!! S*(D^a + D^b)
        pv=matmul(P_v,S)        !!! (D^a - D^b)*S
        F_ab=F_a-F_b            !!! F^a-F^b
        tem=eye-pc              !!! (1-S*(D^a+D^b))
        temp=matmul(matmul(tem,F_ab),pv)          !!! temp=(1-S(D^a+D^b))(F^a-F^b)(D^a-D^b)S
        F_RO=0.5d0*(F_a+F_b+temp+transpose(temp)) !!! F^RO=1/2*(F^a + F^b + temp + temp^T)
        F_ROHF1=matmul(matmul(X_A,F_RO),X_T)      !!! FC=SCe--> F'C'=C'e
        
        eig_vec=F_ROHF1
        call DIAG(eig_vec,eig,NORB)
        C_RO=matmul(X_T,eig_vec)              !!! C=XC'
        call Density1(C_RO,P_a1,NORB,n_el_a)  !!! P_a == D^a
        call Density1(C_RO,P_b1,NORB,n_el_b)  !!! P_b == D^b
        call E_ROHF(NORB,F_a,F_b,P_a1,P_b1,Hcore,E_ele)  !!! E_ROHF == E_UHF obtain energy
        call delta1_P(NORB,P_a,P_a1,delta)    !!!   
        call delta1_P(NORB,P_b,P_b1,delta1)   !!!
        !!! Convergence of judgment 
        if(abs(E_ele-E_ele1)<1.0E-6.and.delta<1.0E-4.and.delta1<1.0E-4)then
        logical_e=.false.
        write(21,*)'****************************'
        write(21,*)'ROHF cycle converged'
        write(21,*)'Number of iterations:',count_iter
        write(21,*)'Electronic energy   :',E_ele
        write(21,*)'Nuclear energy      :',E_nu
        write(21,*)'ROHF Energy         :',E_ele+E_nu
        write(21,*)'***************************'
        write(21,*)    
        call DIAG_1(F_ROHF1,orb_e,NORB)
        write(21,*)"Orbital Energy:      "
        do i=1,NORB
                write(21,'(2X,A3,I2,A1,3X,F16.12)')"orb",i,":",orb_e(i)
        enddo
        !write(*,*) E_ele
        endif
        
        E_ele1=E_ele
        P_a=P_a1
        P_b=P_b1
        !write(*,*)E_ele
        enddo !!!while
        endsubroutine ROHF
        
        
        subroutine E_ROHF(NORB,F_a,F_b,P_a,P_b,Hcore,E_ele)
        integer :: i,j,k,l,NORB
        real*8  :: F_a(NORB,NORB),F_b(NORB,NORB)
        real*8  :: P_a(NORB,NORB),P_b(NORB,NORB)
        real*8  :: Hcore(NORB,NORB),P_t(NORB,NORB)
        real*8  :: E_ele
        P_t=0.0d0
        P_t=P_a+P_b
        E_ele=0.0d0
        do i=1,NORB
           do j=1,NORB
              E_ele=E_ele+0.5d0*(P_t(j,i)*Hcore(i,j)+P_a(j,i)*F_a(i,j)+P_b(j,i)*F_b(i,j))
           enddo
        enddo
        endsubroutine E_ROHF
    
        subroutine Density1(C_new,P_new,NORB,Num)
        !input
        implicit none
        integer :: NORB,Num
        integer :: i,j,k,x
        real*8  :: C_new(NORB,NORB)      !New coefficient matrix
        !output
        real*8 :: P_new(NORB,NORB)       !New density matrix
        P_new=0.0d0
        do i=1,NORB
           do j=1,NORB
              do k=1,Num
                 P_new(i,j)=P_new(i,j)+C_new(i,k)*C_new(j,k)
              enddo
           enddo
        enddo
        endsubroutine Density1

        subroutine delta1_P(NORB,P_old,P_new,delta)
        implicit none

        ! input
        integer  NORB
        real*8  P_old(NORB,NORB)
        real*8  P_new(NORB,NORB)
        integer :: i,j

       ! output
       ! Sum ofmatrix elements square differences
        real*8 :: delta      !Second convergence criterion 


        delta = 0.0d0
        do i = 1, NORB
            do j = 1, NORB
                delta = delta + (P_old(i,j)-P_new(i,j))**2
            enddo
        enddo
        !dsqrt:Square root
        delta = DSQRT(delta / NORB**2)

        endsubroutine delta1_P

